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Abstract. We present and analyse an implicit-explicit timestepping procedure with finite element spa- 
tial approximation for a semilinear reaction-diffusion systems on evolving domains arising from biological 
models, such as Schnakenberg's (1979). We employ a Lagrangian formulation of the model equations which 
permits the error analysis for parabolic equations on a fixed domain but introduces technical difficulties, 
foremost the space-time dependent conductivity and diffusion. We prove optimal-order error estimates in 
the Loo (0, T; L2 (i?)) and L2(0, T; (i?)) norms, and a pointwise stability result. We remark that these 
apply to Eulerian solutions. Details on the implementation of the Lagrangian and the Eulerian scheme are 
provided. We finally report on numerical experiments for an application to pattern formation on evolving 
domains, where we observe novel complex pattern transitions in multicomponent systems and the emergence 
of patterns with different wavelengths under nonuniform evolution. 

1. Introduction 

Since the seminal paper of Turing |29|, time -dependent reaction-diffusion systems (RDSs) have been 
studied as models of pattern formation in natural-process driven morphogenesis and developmental biol- 
ogy (see Murray 1221 for details). An important generalisation of these models consists in considering 
RDSs posed on evolving domains. This stems from the now relatively well-known observation that in 
many cases growth of organisms plays a pivotal role in the emergence of patterns and their evolution 
during growth development |[T3l[22ll . 

RDSs on evolving domains have also a wider scope of application, such as competing- species of 
micro-organisms in environmental biology, chemistry of materials and corrosion processes, the spread of 
pollutants, etc. Numerical simulations of RDSs on growing, or more generally time-evolving, domains 
reproducing empirically observed pattern formation processes are commonly used O [51 El El ED- It 
is essential for scientists to computationally approximate and appreciate the error between simulations 
and exact solutions of such RDSs. Galerkin finite elements ||28]| are among the methods of choice to 
approximate such systems. 

In spite of their widespread use, to the best of our knowledge, no complete error analysis of approxi- 
mating finite element schemes for nonlinear reaction-diffusion systems on evolving domains is available 
in the literature, thus motivating this work. This is a sibling paper to 1 30] where we analysed the well- 
posed nature of (exact) RDSs on evolving domains. In most practical applications, the evolving domain 
is usually a surface embedded in the three-dimensional Euclidean space, but for simplicity we restrict our 
discussion to the case where both the reference domain and the evolving domain are flat, deferring thus 
the analysis of RDSs on evolving curved surfaces. 
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1.1. Problem (RDS on a time-dependent evolving domain). We study a RDS, also considered in (6J 
ITtI . which models a system of chemicals that interact through the reaction terms only and diffuse in the 
domain independently of each other. Given an integer m > 1, the vector u{x^t) G R^, denoting the 
concentration of the chemical species i = 1, . . . , m, at a spatial point x G ftt C R^, d = 1, 2, 3, at time 
t G [0, T], T > 0, satisfies the following initial-boundary value problem 

dtUi{x,t) - DiAui{x,t) ■ [aui] {x,t) = fi {u{x,t)) , xent.te (0,T], 

(1.1) [ly -S/uijix.t) =0, xedQt^tX), 

Ui{x,0) = u^{x), X G 1^0, 



where Qf^ detailed in §2.2[ is a simply connected Lipschitz continuously evolving domain with respect 
to t G [0, T], and D := (I^i, . . . , Dm)^ is a vector of strictly positive diffusion c oeffi cients. Detailed 



assumptions on the nonlinear reaction vector field / := (/i, . . . , /^)^ are given in ^2.5 The convection 
a = (ai, . . . , adY is induced by the material deformation due to the evolution of the domain. The initial 
data is a positive-entry bounded field. Since we are primarily interested in pattern formation phenom- 
ena that arise as a result of self-organisation within a domain without outside- world communication we 
consider homogeneous Neumann boundary conditions, but other types of boundary conditions could be 
studied as well within our framework. 

1.2. Main results. The core result in this paper is Theorem |5 . 1 1 where we prove optimal convergence 
rates of the discrete solution in Loo(0, T; L2(r^))^ and L2(0, T; H^(r^)^) (where 1] is a transformed 
version of Vtt to be described next). Our theoretical results are illustrated by numerical experiments, aimed 
mainly at quantifying the pattern formation phenomena related to the type of growth in the domains. The 
results of numerical experiments exhibit the emergence of complex patterns when multispecies RDSs are 
considered as well as patterns of different wavelengths when nonuniform domain evolution is considered. 
To the best of our knowledge, this article is the first to expose such results. 

1.3. A Lagrangian approach. We employ, both for the analysis and the implementation of the compu- 
tational method, a Lagrangian formulation of Problem ] 1.1 [ in the sense employed in fluid-dynamics, i.e., 
where the evolving domain, Q.t ^ R^/i^ the image of a time-dependent family of diffeomorphisms A^t on 
a reference domain Vt ^ R^. The m parabolic equations with constant diffusion coefficient constituting 
the RDS on Vtt are thus pulled-back into equations on a fixed domain, albeit with space-time dependent 
coefficients. The fixed domain setting permits us to use the standard Bochner space machinery needed for 
evolution equations of parabolic type. On the other hand, we are thus left to deal (computationally and 
analytically) with three interacting difficulties: (1) m coupled equations, (2) the nature of the nonlinearity 
/ coupling the equations, (3) the non-constant diffusion and velocity coefficients, especially as functions 
of time. Our approach in tackling the nonlinearity consists in constructing a suitable globally Lipschitz 
extension of the nonlinear reaction field that coincides with it in a neighbourhood of the exact solution 
and then proving that both the exact solution and the numerical solution are confined to the domain of the 
original (non-extended) nonlinearity. We use mainly parabolic energy techniques, but must have some 
pointwise control in order to bound the nonlinearities. 

Although all our error estimates are derived for the Lagrangian formulation, given that the domain 
evolution is prescribed, they carry in a straightforward manner to the Eulerian framework. The situation 
would be more delicate if the domain evolution was itself an unknown as a geometric motion coupled 
with the RDS, but this is outside the scope of this study. 

1.4. Implicit-explicit schemes. The fully discrete method that we analyse, is a fully practical method, 
implemented in the ALBERTA toolbox (code available upon request), using an implicit-explicit back- 
ward Euler scheme to derive the time-discretization 1 18|. On fixed domains, Zhang et al. L32J analyse a 
second order implicit-explicit finite element scheme for the Gray-Scott model and Garvie and Trenchea 
191 analyse a first order scheme for an RDS that models predator prey dynamics. Recently Mackenzie 
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and Madzvamuse 1 16] analysed a finite difference scheme approximating the solution of a linear RDS 
on a domain with continuous spatially linear isotropic evolution. Our study is novel, in that, we pro- 
pose and analyse a finite element method to approximate RDSs on a domain with continuous (possibly 
nonlinear) evolution. This creates a space-time-dependent coefficients impacting the diffusion and the 
time-derivative term which complicates the fully discrete scheme's analysis and requires a careful treat- 
ment of the timestep, depending on the rate of domain evolution. In spite of it being only first order in 
time, the proposed implicit-explicit method is robust for the applications we have in mind, where long 
time integration is essential and the problems are often posed on complex geometries such as the surface 
of an organism. 

1 .5. Outline. The structure of this paper is as follows: in ^we introduce the notation employed through- 
out this article, we state our model problem together with the assumptions that we make on the problem 
data and the domain evolution. We present the weak formulation of the continuous problem and de- 
fine a modified nonlinear reaction function which we introduce for the analysis. In ^ we present the 
semidiscrete (space-discrete) and the fully discrete finite element schemes with some remarks regarding 
implementation, allowing the practically minded reader to skip over the analysis through to ^ We then 
analyse the semidiscrete scheme in Q and the fully discrete scheme in ^ proving optimal rate error 
bounds as well as a maximum-norm stability result, whereby the stabilising effect of domain growth ob- 
served in the continuous case is preserved at the discrete level and in the numerical schemes. In ^we 
provide two concrete implementations (Eulerian and Lagrangian) of the finite element scheme with a set 
of reaction kinetics commonly encountered in developmental biology, considering domains with spatially 
linear and nonlinear evolution. In ^ we present computational experiments to illustrate our theoretical 
results and to study as of yet unexplored solution behaviour exhibited by RDSs on evolving domains in 
two- space dimensions. 



2. Notation and Setup 



In this section we define most of the basic notation for the rest of the paper, introduce the evolving 
domain framework, set the detailed blanket assumptions and introduce a pulled-back version of Problem 

o 

2.1. Calculus and function spaces. Given an open and bounded stationary domain IT c [R^ and a 
function 77 G C^(II;[R^),we define the Jacobian matrix of 77 



(2.1) 



Vr/(x) := 



, for X G n, 



and the divergence of 77 
(2.2) 



V • ri{x) := ^dx^r]i{x). 



For the case of a scalar- valued function r] G C^(II; R), we define the Laplacian of r] 



(2.3) 



Ar]{x) := "^d^.^.rjix) 



i=l 



In an effort to compress notation for spatial derivatives, we introduce the convention used above, that 
if the variable with respect to which we differentiate is omitted, it should be understood as the spatial 
argument of the function. 
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We denote by Lp (H), W^'^ (E) and (E) the Lebesgue, Sobolev and Hilbert spaces respectively, 
as defined by the following: for measurable rj and for p,k& [1, oo) 



(2.4) Lp (n) :={r]: 



(2.5) Loo (n) sup \r]{x)\ < oo 



(2.6) w'^'P (n) := i 7? e Lp (n) : ^ v"7? e Lp (n) i , 

[ a</c J 

(2.7) H^(n) := W^'^(n). 

For measurable functions 7^, /i : 11 ^ [R, we introduce the following notation 
(2.8) 

Ju 

(2.9) 11-11- • • — 





:= / r]{i 




Ju 


(n) 


•= iv^v)} 




:= ||VS 


^IIh^ (n) 


:= (ll^ll: 



1/2 



(2.10) 

(2.11) ll^llHMn):= \^\v\\l(u)+J2\v\w(n) 
For vector valued functions rj^ fi .U ^ R^, we denote 

(2.12) (^,M)n-'=y^/ 

with the corresponding modifications to the norms and seminorms. 

2.2. Evolving domain. Let ^7 C be a simply connected domain with Lipschitz boundary; we will 
call it the reference domain. We define the evolving domain as a time-parametrised family of domains 

(2.13) {Vtt '= ^t{^)}o<t<T where : ^ is a C^-diffeomorphism for each fixed t G [0, T]. 
The Jacobian matrix of v4t(-), its determinant and inverse will be respectively denoted by 

(2.14) J,(0 := VA(0,^t(0 '= det J,(0 andK,(0 := [VA(0]"' breach (C,t) e^x [0,T]. 
We will use also the evolution induced convection on the evolving domain 

(2.15) a{x,t) := dtAt{A[^{x)) forx G Qt and t G [0,T]. 
From classical results 111 we have the following expression 

(2.16) dtJ (C, t) = J,(0 V • a for (C, t) e Cl x [0, T]. 

To aid the exposition we define Q to be the topologically cylindrical space-time domain: 

(2.17) Q:={{x,t):xent.te[0,T]}. 

We now introduce notation to relate functions defined on the evolving domain to functions defined on the 
reference domain. Given a function g : Q ^ R we denote by ^ : Cl x [0, T] [R its pullback on the 
reference domain, defined by the following relationship 

(2.18) 9{^,t):=g{At{^),t) i^,t) € Cl x [0,T]. 
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Assuming sufficient smoothness on the function g, using (2.18) and the chain rule we may relate time- 
differentiation on the reference and evolving domains 

(2.19) dtgi^.t) =829 + ' V^] , {^.t)enx [0,T]. 

The right hand side of p.l9| ) is commonly known as the material derivative of g with respect to the 
velocity a. The following result relates the norm of a function g : Q ^ [R on the evolving domain with 
its pullback g on the reference domain: 

(2.20) \\9\\lin,) = {Jtg,g)n=-\\9\\l- 
For the gradient of a sufficiently smooth function ^ : Q ^ [R, we have 

(2.21) W'^gWlin.) = {JtKtVg,KtVg)^ = {BtVg,Vg)^ =: \g{;t)\l^, 
where B := JKK^ . For t e [0, T] we define the bilinear form 

(2.22) ht {v, w) := {BtVv, Vu;)^ for v,w e 



Assumption 2.3 implies that there exists /i, /i G [R+ such that for z = 1, . . . , m and all v eH^ (Cl), 



(2.23) ^^\\^i\\l^ci)<btii,v) < l|S||L^(f,)||Vt)||^^(^) =M||Vt)||2^(^). 

2.3. Assumption (Regularity of the mapping). We assume the following regularity on the mapping A. 
between the the reference domain Cl and the evolving domain Qt • 

(2.24) Ae&{nx[0,T]) and At e C^^^ (Q) for t e [0,T], 

where k will be taken equal to the degree of the basis functions of the finite element space defined in the 
following section. To ensure the mapping is invertible we assume the determinant of the Jacobian J of 
the mapping A (cf. ( 2.14 )) satisfies 

(2.25) J >Omnx[0,T]. 



2.4. The RDS reformulated on the reference domain. Using ( |2.18[ ) — ( |2.21[ ) and a change of variable 
in the divergence, we obtain the following equivalent formulation of Problem 
Denote by u : ft x [0, T] the function that satisfies for z = 1, . . . , m, 



1.1 on a reference domain. 



dtUi-^W ■{B\/ui)^Ui\/ -0. = fi{u), onClx (0,T], 



(2.26) 



i> ■ BVui = 0, ondVt x (0,T], 
where 

(2.27) a(C,t) = dtAt{$,) for (C,t) G x [0,T]. 

2.5. Assumption (Nonlinear reaction vector field). We assume throughout that / is of the form 

(2.28) /i(^) = for all z G Dom f =' I and each z = 1, . . . , m, 
for some vector field F G C^(/). As a result / G C^(/) and it is locally Lipschitz. 

2.6. Ass umption (Existence and regularity). We assume the global existence of a solution u to Prob- 
Furthermore we assume u is in H^~^^{tt)'^ with dtu in H^~^^{tt)'^ where £ is the polynomial 



2.4 



lem 

degree of the finite element space defined in the following section. 



When confusion may arise, such as in j2.19t ) for a positive integer i, dif denotes the partial derivative with respect to the i-th 
argument of the function /. When there is no risk of confusion we write dtf for the time derivative of a time-dependent function / 
even when such a variable is not explicitly written in the arguments. 
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2.7. Remark (Applicability of Assumption |2.6| l. In |30|, we proved the global existence of positive 
classical solutions to Problem [LT] for a class of RDSs with positive initial data on domains with bounded 
spatially linear isotropic evolution. In this case, if / belongs to C^(R!p), the solution u is shown to 
be contained in a rectangle Yli=i[u^ ^u^]^ ^i^^ u~ > and assumption 2.5 holds with the region 



h •= nti K - + ^] a subset of 

2.8. Weak formulation. For the purposes o f con structing a finite element discretisation, we introduce 

The problem is to find Ui e L2 (O, T; {Cl)) with 

D,bt {S/u^,\/x) = {Jf^{u),x)c^^ 



a weak formulation associated with Problem 

dtUi e L2 (0,T;H"^ {Cl)) such that 

(2.29) / J (^dtUi ^UiV-a {At{^),t 



2.4 



for all X ^ H {Q). Using the expression for the time-derivative of the determinant of the Jacobian (2.16), 
we have 

(2.30) {dt{Jui),x)^+Dibt{Vui,Vx) = {JMu),x)ci, Vx e (fi). 

We shall use ( |2.3Q| ) to construct a finite element scheme to approximate the solution to Problem [LT] on 
the reference domain. In ^we illustrate that the resulting scheme may be solved both on the reference 
(Lagrangian scheme) or the evolving (Eulerian scheme) domain. 



2.9. Extended nonlinear reaction function. In general the techniques used to show Assumptions 2.5 



and |2.6| holds utilise the maximum principle |22l EOI • In the discrete case, the absence of a maximum 
principle | 28, pg. 83] means we cannot guarantee the discrete solution remains in the region !§. For the 
purposes of our analysis we intro duce a modified globally Lipschitz nonlinear reaction function. Recalling 
Is and F from Assumption 2.5 we define F e C^(W^) such that 

F(z) 
f\z) 



(2.31) 



= F{z), forz€h 
<C, for^GR™, 



and fi{z) := ZiFi{z), for z G R*^ 



The function F is guaranteed to exist due to Assumptions 



2.5 



2.6 



Th. 1, §6.5] and the use of an appropriate cut-off factor. If it is a solution of ( |1.1[ > 
(2.32) f{u) = f{u). 

Thus, we may without restriction replace / with / in \\.\\ . 

3. Finite element method 



the Whitney Extension Theorem (H 



In this section we design the finite element method, first by discretising Problem 2.4 in space only, 
discussing some properties of the semidiscrete scheme and then passing to the fully discrete scheme. 

3.1. Spatial discretisation set-up. We shall split the spatial and temporal discretisation of Problem [TTT] 
into separate steps. For the spatial approximation, we employ a conforming finite element method. To 
this end, we define ^ a triangulation of the reference domain. We shall consistently denote by /i := 
max^^^ diam(5) the mesh-size of We assume the triangulation ^ fulfils the following properties: 

• By 5 G ^ we mean s is an open simplex (interval, triangle or tetrahedron for d = 1, 2 or 3 respec- 
tively). 

• ^ is conforming, i.e., for any s^k e sDkis either 0, a vertex, an edge or a face common to s and 
k or the full simplex s = k. 



No error due to boundary approximation, i.e., U 



Cl (we make this assumption for ease of 



exposition and it may be relaxed depending on the application). 



Implicit-explicit timestepping with finite element approximation of reaction-diffusion systems on evolving domains 



7 



For a sequence {^}^^ of conforming triangulations, we assume the quasi-uniformity of the sequence 
holds, i.e., there exist Ci, C2 independent of i such that 



(3.1) 



Cih <hs< C2K, for alls G 5^, z = 1, 2, . . . , 



where /i° and hs are the radius of the largest ball contained in s and the diameter of s respectively. 
Furthermore we note that our assumption of quasi-uniformity implies that the family of triangulations is 
shape-regular [26, pg. 159]. 

Given the triangulation we now define a finite element space on the reference configuration: 



(3.2) 



:=|l>Gi^^(r^):l>|sis piecewise polynomial of degree ^| . 



We utilise the following known results about the accuracy of the finite element space V. By the definition 
of V, we have for v G H^^^{Q) (see for example Brenner and Scott |4| or Thomee |28 |), 



(3.3) 



inf < \\v 
l>ev 



^1 



L2 {n) 



h\\V{v-^) 



IL2 {n) 



} < Ch'^^ 



In the analysis we shall make use of the fact that ( |3.3| ) is satisfied by taking the Lagrange interpolator 
: H^+^ (Cl) ^ V in place of 4>. Let : C° ^ V be a Clement type interpolant 15J and ^ + 1 > f , 
where d is the spatial dimension. The following bound holds 



(3.4) 



We shall make use of the following inverse estimate valid on quasiuniform sequences of triangulations: 



(3.5) 



l^ll 




3.2. Semi discr ete approximation. We define the spatially semidiscrete approximation of the solution 
to be a function : [0, T] V, such that for z = 1, . . . , m, 

dt{Ju^), 4>)^ + (d.bvu'i, v4>)^ = ( J^(n^), 4>)^ , vl> G V, 

where A^ is the Lagrange interpolant. 

3.3. Proposition (Solvability of the semidiscrete scheme). Let Assumptions \2^ and \2.3\ hold. Then, 
the semidi scre te scheme \ 3. 6) possesses a unique solution G L^o (0, T)^ 

Proof. In 

equations 

and the construction of / ( |2.31| ), we have that J, / and their product are continuous globally Lipschitz 
functions. From ODE theory (for example 11241 ) we conclude that ( |3.6| ) possesses a unique bounded 
solution. □ 



3.6) if we write Ui{t) as ^^^^ ^j^j^ we obtain a system of dim(V) ordinary differential 



or each i. By assumption the initial data for each ODE is bounded. From Assumption 2.3 



3.4. The effect of domain evolution on the semidiscrete solution. We now examine the stability of 
( |3.6| ) and show that domain growth has a diluting or stabilising effect on the semidiscrete solution, mir- 
roring results for the continuous problem 1 14|. Taking ^ = in (3.6) gives for z = 1, . . . , m, 

(3.7) (dtiJul), ul)^ + Dih (Vut Will) = (jfiiu"), ill 



For the first term on the left of (|3.7|) we have 
(3.8) 
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Application of Reynold's transport theorem fl] gives 



(3.9) 



Ci 2 \dt 



Now dealing with t he right hand side of (3/7(, using ( 2.31| l and the mean- value theorem (MVT) we have 
with C from t..3l\ 



(3.10) 

Therefore we have 
(3.11) 



< 



/i(0)|+C^|A,^| 



<c{jj2\ui\, 



Applying Young's inequality gives 
(3.12) 

where (o) ^ "^^ depends on /i(0) 



m + 1 , 



■^/i(0)' 



(3.13) 



E 



. Summing over i we have 
1 



< [Cm 



- h\\2 



Using ( [22T] ), ( [3^ and ( [313] ) in ([TT]) gives 

d, 



dt 



\v!'fj^ + 2 ^ Al^f 1^ < ( J [2Cm + 1 - V • a 



(3.14) 

Finally, integrating in time and applying Gronwall's lemma we have 

(3.15) < (||n^(0)|p^.+2tCy(o))exp sup |2Cm + 1 - V • a ^) | M • 

^ ^ \j7x[0,T] ^ ^ / 

From (2.16), the dilution term V • a has the same sign as dtJ and is therefore positive (or negative) if 
the domain is growing (or contracting). Thus, domain growth has a diluting effect on the L2 (^^t)^ norm 
(c.f., ( [2.2Q| )) of the solution. 

3.5. Fully discrete scheme. We divide the time interval [0,T] into N subintervals, = to < ... < 
In = T and denote by := tn — tn-i the (possibly nonuniform) time step and r = max^ r^. We 
consistently use the following shorthand for a function of time: := f{tn), we denote by df^ := 

For the approximation in time we use a modified implicit Euler method where linear reaction terms and 
the diffusive term are treated implicitly while the nonlinear reaction terms are treated semi-implicitly us- 
ing values from the previous timestep (a 1-step Picard linearisation). Our choice of timestepping scheme 
stems from the numerical investigation conducted by Madzvamuse 1 181. 
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The fully discrete scheme we employ to approximate the solution of Problem 1.1 is thus, find U^^ G V, 
forn = 1, . . . , A^, such that for z = 1 . . . , m, we have 



(3.16) 




DJ[BVU,r,V^ 



C/° = A 



where A'' is the Lagrange interpolant and is as defined in (2.28 (. 

3.6. Moving mesh formulation. Alternatively, and in a more physically intuitive way, we may look to 
approximate the solution to on a conforming subspace of the evolving domain. To t his e nd we define 
a family of finite dimensional spaces V", n= [0, . . . , A'^] such that with $ as defined in (3.2): 



(3.17) 



which also defines the triangulation ^",n = [0 , . . . , A''] on the evolving domain. Using (3.16) and (3.17) 
we have the following equivalent finite element formulation on the evolving domain: find (7" G V", for 
n = 1, . . . , A^, such that for i = 1 . . . , m. 



(3.18) 



d 



-A {VC/f,V*")f,^„ 

where A'' : ($7o) -> Vq is the Lagrange interpolant. By |2.19| ) and |3.17 ) 
(3.19) dtHO = [dt^ + a • V$] {At {$,)) = 0. 



4. Analysis of the semidiscrete scheme 

We now prove that the semidiscrete solution converges to the exact one with optimal order in the 
Loo (0, T; L2 {n)"^) norm and the L2 (0, T; ((^)^) seminorm. 

4.1. A time-dependent Ritz projection. A central role in the analysis is played by the Ritz, or elliptic, 
projector, defined, as in Wheeler (STJ, for each t G [0,T], by • H^{^) V such that for each 



(4.1) 
(4.2) 



6, (i),4>) =bt(^Rtv,^) V4>G V, 

and / [Rtv - v] = 0. 
Jq 



The constraint (4.2) ensures Rt is well defined. Differentiation in time in (4.1 ) with v = Ui yields 

(4.3) b (dt{ui - Rui), $) + {{dtB)V{ui - Rui),V^)^ = 0, V$ e V. 

To obtain optimal error estimates, we now decompose the error into an elliptic error (the error between the 
Ritz projection and the exact solution) and a parabolic error (the error between the semidiscrete solution 
and the Ritz projection): 



(4.4) 



u= {u^ - Ru) + [Rtii - u) 



where the equality defines p = (pf^, . . . , p^)^ and s = (^i, . . . , £m)^- 
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4.2. Lemma (Ritz projection error estimate). Suppose assumptions \2.^ and \2.3\ (with k = £) hold and 
let R be the Ritz projection defined in (4.i] ). Then the following estimates hold: 



(4.5) sup <! \\RMt) - iCi)r^ + E 11^ (^^^^W - ^^W) 



te[o,T] 



1=1 



< C{A,u)h 



2(^+1) 



(4.6) 



sup \\dt {RMt) - m) \\l (^)^ + E 11^^^ i^Mt) - Mt)) \\l (^) < C{A, u)h 



^2(^+1) 



tG[0,T] 



i=l 



Proof. Using ( |2.23 ) and \AA) we have for i = 1, . . . , m 
(4.7) 



VIIL2 (^7) 



Il2 (^7) P*Ih^+1(^^) ' 



which shows the energy norm bound of (4.5 ). To show the L2 estimate we use duahty. Fix at G (0, T] 
and consider the solution ip of following elliptic problem 



(4.8) 

Note that 



V • {Bt\/2p) = 4) in ft, BtV^ 'U = {)ond(l, 2p = 0. 



< C 



Hi(^) 



L2 m 

(4.9) 

We therefore have 
(4.10) 

Furthermore we have the estimate 



as for any v 



L2 m 



<C\vU 



Hi(^) • 



(4.11) 



H2(n) 



< C\m\^^ < C\\B A^W^^ ^^,=CU^\/-B- V^IIl, < C 



IL2 (^7) 



L2 (^^) 



IL2 



Here we have introduced the notation that the divergence of the tensor S is a vector defined such that for 
z = 1, . . . , (i 



(4.12) 



i=i 



Thus testing (4.8 ) with ii and using (4.1 ( we have 

(4.13) \ /n V J 

< /S||V£i||L,(o)l|V(^- J'^^)||L,(n) < Ch'^^' l«lH^+Mn) l^lH=(n) < Ch'+\ 



which completes the proof of ( |4.5| ). For the proof of ( |4.6| ) we have using ( |4.3| ) and the fact that the gradient 
commutes with the time derivative (as we work on the reference domain) that for z = 1, . . . , m, and for 
each 4> G V, 



(4.14) 



I^W^dtiiWl^^^^ < bt {dtii.dtii) = h {dtSi,^ - dtu^ + h (dtii^dtRm - ^ 
= bt (^dtii, 4> - dtu?j + (dtBVsi, V(4> - dtRtUi) 
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Taking 4> = XhdtUi in (4.14) gives 



(4.15) 



llL.(Q) + l|v(x''5tii,-5tU,)| 



L2(n) 



IL2 (^^) 



where we have use Young's inequaUty in the final step. The previous estimate \A.5) completes the proof 
of the energy norm boun d in \A.6\ . For the L2 estimate we once again use duality. Testing problem ( |4.8| ) 
with dtii and using (4.3 ), we have for z = 1, . . . , m, and any <l> G V 



(4.16) 



Taking $ = XhdtUi in (4.16) gives 













<c 




H2(n) I 







ft./Li||V9t£ 



(4.17) 



/i||5t-B| 



Loo (n)ll^'^«llL2 (O) 



|5tS|lL..olkillL.(Q) 



l|V^|lL,(ml|V5tB||L„ 



where we have used integration by parts to estimate the last term in (4.16). The previous estimates and 
Assumption |2. 3 1 complete the proof. □ 



4.3. Theorem (A priori estimat e for the semidiscrete scheme). Suppose Assumptions 2.5 and\Z 



hold. Furthermore, let Assumption 2.3 hold (with k = i). Finally let u be the solution to Problem Ua 



Then, the following optimal a priori error estimate holds for the error in the semidiscrete scheme: 
(4.18) 



sup {\\u''{t)-umlrCi)^} + Y. 

t€[0,T] ^ '^'^"^ J ^ Jo 



Proof. Using the decomposition ( |4.4| ) and Lemma 4.2 we hav e a bo und on the elliptic error and it simply 
remains to estimate the parabolic error p^. To this end, we use ( 3.6 ) to construct a PDE for by inserting 



in place of and taking ^ = p^. Using (2.21 ) we obtain for z = 1, . . . , m 



(4.19) U (Jp'l) +A|Vpf||j = {^h{u''),Jp'})^-(dt{JRtUi),p'})^-bt ivRtUuVpl) . 



Using (2.30(, (2.32( and (4.1 1 gives 



(4.20) {dt (Jp^) , p'lj^ + Di\Vp1\l = { Mu'') - Mu),Jp1 



dt{Jei),p'} 



Dealing with the first term on the left of ( 4.20| l as in ( |3.9[ >: 

_ 1 / d 



(4.21) 



dt{jp'l),p'l 



2 Vdt'"^''"'+\'^'^''^'''^'^*^^^^^'^^o 
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Dealing with the first term on the right of ( 4.20| l using ( |4.4| i and the MVT we have 



(4.22) 



.J 



Pi 



Applying Young's inequality: 



(4.23) 



Summing over i we have 



/ n 



<C 



(4.24) 



E 



<c\ 



Dealing with the second term on the right of (4.20): 



(4.25) 



dt {Jei),p\ 



< 



< 



where we have used Young's inequality for the second step. Now using (2.16) and summing over i we 
have 



E 



(4.26) 



< 



1 



Jph V ■a{At{£,),t) 



WdtiVj^ + \\dtJ\\ 



if 



Loo (nx[o,r])ll'^llL2 (Q) 



Combining ( |43T] l, ^2^ , 



(4.27) 



dt 



IP IIj" 



::||2 



where we have used the fact that Assumption 2.3 implies J^dfJ G Loo x [0^ ^]) • Integrating in time, 

onwall's L( 

Ea/ 

i=l -^0 



using Lemma [43| and applying Gronwall's Lemma we have 



(4.28) 



lirWII 



|V/5; 



/i|2 

i Is 



< c 



(l|p^0)||2,„+/.^(^+i)). 



To estimate p^(0), we note 



(4.29) 



< ||n(0) - A^n(0)|| + ||s^|| < Ch'+\ 



where we have used ( |3.3|) , the assumption on the regularity of the exact solution and Lemma [4^ in the 
last step. Assumption |2. 3 1 and the equivalence of norms p.2Q| ) completes the proof. □ 

5. Error analysis of the fully discrete approximation 



In this section we provide the convergence result for the fully discrete scheme ( 3.16 ). The main result 
of this paper is Theorem |5.1| whose proof with give in details. We follow that up with a convergence 
result in the Loo (^) norm which allows the use of the original / (without extending to / in the numerical 
method). 
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5.1. Theorem (A priori estimate for the f ully discrete scheme) 

(with k 



2.3 



Suppose Assumptio ns \2.5\ and 
£) holds. Let U be the solution to (3.i6|). Finally, 



hold. Furthermore, suppose Assumption 

suppose the timestep satisfies a stability condition defined in (5.ii| ). Then, the following optimal a priori 
estimate holds for the error in the fully discrete scheme: 



(5.1) 



<c(a^,5) + r') , forne [0,...,Ar], 



5.2. Remark (Error estimate for the evolving domain scheme). The schemes ( |3.16| ) and ( |3.18| ) are 

equivalent up-to numerical quadrature. Thus Theorem 5J_ also provides an error estimate for the evolving 
domain based scheme ( |3.18| ). 

Proof of Theorem |5.1 Decomposing the error as in ([4.4|) we have 



(5.2) 



L2 {^Y 



From Lemma 4.2 we have the following bound on the elliptic error: 



(5.3) 



||2 



IL2 {^lY 



< forn G [0,...,A^]. 



Therefore it only remains to estimate . Constructing an expression for as in (4.19 ), using ( 3.16 ) and 
\A.\) we obtain for z = 1, . . . , m, 



(5.4) 



{d[ji,r, p^)^ + ( (a - dt) [ju^r, p?) , , 



where we have used (2.30) for the second step and F is as defined in (2.31 ). Using Young's inequality 
for the first term on the left hand side of (|5.4|) gives 



(5.5) {d[Jpir,P7)ci> 



where we have used (2.20). Summing over i we have 



(5.6) 



n-l 



i=l 



2" J" 



2t„ 
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Using 5.2 and the MVT for the first term on the right of (|5.4|) gives 



(5.7) 



-n-1 



Tndu] 



n-1 1|2 



L2 (^^)" 



where we have used Young's inequaUty for the second step. Summing over i we have 



E 

2=1 



n-1 



(5.8) 



jn 



IIJ^ 



l|2 



-n-1 1|2 



Applying Young's inequality to the second and third term on the right of (|5.4|l gives 



(5.9) 



(d-dt) [Juir,pi 



< 



o7fj 



Using ( [5^ , and in pA} gives 



1 



1 



(5.10) 



< 



2 II Jn llLoo(^^) 

1 _ 7^ 

+ ^^ifeiii 



B 



2t„ ' --"jn-i^L^m 

c5||J"|lL^(n)(ll^"ll 

1„ 1 „ 



|p"-'ll5" 



2 

L2 (^^)'^ 



-n-1 1|2 



IL2 (n)- 



l-"^«"llL(n). 



Let > be such that, for r < r' and forn = 1, . . . , A^, 



11L2 (^7)^ 



0^ 



(5.11) 
Such a r' exists since 

(5.12) 



1 



lim 



2" Jn llL^(n) 



CCt > 0. 



1„ J" 



+ CCr 
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For r < r\wQ have 

m 

(5.13) WpXj^ + ^CrAIVpnlj < C ((7"||p"-i% + tTI") , 



where C" = 1 + TC\\y{^ 
(5.14) 



and 



R-:=C||J''||,_|S,UI«iiL.<0,. 



-n-1 1|2 



11L2 {nr 



\\rdu 



n\\2 



+2llj;rllL.(o)(llWI 

Therefore, forn = 1, . . . , A/", 



L (0)^ + 11 



L2 (^^) 
L2 (^^) 



(5.15) 



i=l 



<c{ llc'^\\prj^ + rJ2ll^'n^)- 

j=l i=j / 



. /c=l 



Considering the first two terms on the right of (5.14), we have forn = 1, . . . , A/" 



c\\r 



xn-1 1|2 



11L2 (n)-' 



(5.16) 



1^^112 



<2C sup ||J'||l^(^)||c II 



where we have used Assumption 2.3 and Lemma 4.2 DeaHng with the third term on the right of ( 5.14 ), 
we have 



(5.17) 



c\\r\ 



Loo i^) 



\rdu 



n\\2 



c\\n 



L2 (Q)^' 



where we have used Assumptions 2.6 and 2.3 For the fourth term on the right of ( 5.14 ) we have 



2 II jn IIloo m 



\\d[Js] 



< 



1„ 1 



5t[Je]Ms| 



L2 



(5.18) 



< C sup 



IL2 (^7)- 



where we have used Assumption 2.3 for the second step and Lemma 4.2 for the final step. Finally, for the 
fifth term on the right of ( |5.14| ) we have 



'J 



nllL..(^)ll(^-St) [JU] 



n\\2 



IL2 (O)-^ 



(5.19) 



- II jn IlLoo i^)-^^ 

< Ct'^ sup [Wdtii 



-/* (s-e-^)dtt[JuYM 



2 

L2 (Ci)- 



s\\2 

L2 iCiY' 



s\\2 



IL2 (nr 



16 



O. LAKKIS, A. MADZVAMUSE, AND CHANDRASEKHAR VENKATARAMAN 



where we have used Assumption 2.3 for the second step and Assumption 2.6 for the final step. Combining 
15161 ), ( [5T7l ), ( [5l8l ) and ^J9^ we have 

(5.20) < C (h^^^^^^ + r^) forn = 1, . . . , AT. 
Using ( |4.29| ) we have 

(5.21) ||p°||5„ = 11/(0)113™ <C/iW). 

Applying estimates ( |5.20| ) and ( |5.21| ) in ( 5.13| ) completes the proof of Theorem [51] □ 

5.3. Remark (Stability of the fully discrete scheme). The timestep restriction ( |5.11| ) is composed of a 
term arising from domain growth (the term involving the determinant J of the diffeomorphism A) and 
a term arising from the nonlinear reaction kinetics (the term containing C). It is worth noting that for 
a given set of reaction kinetics, i.e., a given C, larger timesteps are admissible on growing domains (as 
we have ||-^^~^/-^^ ||l (Ci) ^ consider for illustrative purposes the heat equation, i.e, the case 

C = 0, we recover unconditional stability on growing domains whereas for contracting domains (5.11 
implies a stability condition on the timestep dependent on the growth rate. 



5.4. Remark (Qualitative estimates on the exact solution). In practise only qualitative a priori esti- 
mates are generally available for the exact solution and the region Is defined in Assumption |2.6|is not 



explicitly known. Thus, we cannot construct the function / defined in (2.31 ). To this end, we introduce 
the following assumption to circumvent the construction of /. 



5.5. Assumption (Dimension dependent polynomial degree). We wish to invoke estimate ( |3.4[ ) with a 
positive power of h and thus we require the degree of the finite element space to satisfy 

d 



(5.22) 



£> 



1, 



where d is the spatial dimension. Thus, we require piecewise linear or higher basis functions for d < 2 
and at least piecewise quadratics for d = 3. 



5.6. Remark (Maximum-norm bound of the discrete solution). Let assumptions in |5.5| and Theorem 
I5.1lbe valid. Then we have 



(5.23) 



U 



and for sufficiently small mesh-size h the discr ete so lution U to Problem ( 3.16 ) is in the region Is for 
all n G [0, . . . , A^]. Thus, we may replace F in (3.16) by F. 

Indeed, forn G [0, . . . , A'] we have for the Clement interpolant 

(5.24) \\vr-U"h .o^^ < ||X^n"-Z7"||, ,o^^ + - X^u^ 



Using ( 3.4 ) and ( 3.5 ) gives 



(5.25) 



Error bound (5.23) now follows from (3.3) and Theorem 5.1 Thus, if h is taken sufficiently small we 
have 

5 
2' 



(5.26) 



sup \\u'' -U IIl^(^)^ < 



ne[0,...,N] 
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Therefore, U e I5 for all n G [0, . . . , N] and thus, f{U) = f{U). 



6. Implementation 

In this section we illustrate the implementation of the finite element scheme with explicit nonlinear 
reaction functions. We demonstrate that scheme ( |3.16| ) on the reference domain or scheme ( |3.18| ) on the 
evolving domain give rise to equivalent linear systems. To illustrate a concrete application of the proposed 
scheme we consider the following widely studied set of reaction kinetics. 

6.1. Definition (Sciinakenberg's "activator-depleted substrate" model |[T0l[l5l|25|). We consider the 
following activator depleted substrate model, also known as the Brusselator model in nondimensional 
form: 



(6.1) 



where < a, 6, 7 < oc. 



In matrix vector form scheme ( |3.16| ) equipped with kinetics ( |6.1| ) and appropriate initial approxima- 
tions W2 is: To solve for W^^ ^2 ' ^ = [1, . . . , A^], the linear systems given by 



(6.2) 



+DiS +jN.]W. =^M W. 



-M 



D2S +-/N2)W2 



-^M Wn 



■jbF , 



where Wi and W2 represent the nodal values of the discrete solutions corresponding to iii and 112 re- 
spectively and the equations are nondimensional such that either Di or D2 is equal to 1 . The components 
of the weighted mass matrix M, the weighted stiffness matrix S and the load vector F on the reference 
frame are given by 

(6.3) 

(6.4) 
and 
(6.5) 



:= /_ J"*„$;3, 



on 

^a(3 - 



For reaction kinetics (6.1 ) the components of the matrices arising from the Picard linearisation A^i and 
N2 are given by 



(6.6) 

and 

(6.7) 



dim(V) dim(V) 



dim(V) dim(V) 



respectively. 

We now illustrate the implementation of scheme ( 3.18 ) where the linear systems are assembled on the 
evolving domain. By the definition of V"^ (3.17 ) we obtain the following time dependent mass {M'^) and 
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Stiffness (5^) matrices if we assemble the linear systems on the evolving domain 

(6.8) := f $»$^ = M^p, 

(6.9) S^p := f • = S2p, 
and the load vector on the evolving domain is given by 



(6.10) 



n pn 

a a . 



We thus obtain the following linear systems: 



(6.11) 



where for reaction kinetics ( |6.1| ) the components of the matrices arising from the Picard linearisation 
Ni = Ni are given by 



(6.12) 



dim(V) dim(V) 
r]=l 'd=l -^^t^ 

with analogous modifications for N2 • 



Both formulations (6.2) and (3.18 ) result in the same linear algebra problem. Solve for vectors 6^, i 



(6.13) 



A"6^ = c^-\ forn = l,...,iV. 



The matrix is symmetric sparse and positive definite. We therefore use the conjugate gradient (CG) 
algorithm 1 12] to compute the solution to the linear systems. 




Figure 1 . A simple example of the reference and evolving domain with the associated 
mapping, mesh-size and triangulations. 
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6.2. Remark (Quadrature on the evolving domain). As we do not have to compute the Jacobian of 
the mapping, assemblage of the Hnear systems is faster on the evolving domain. However, in the previous 
analysis we have neglected errors due to variational crimes, such as the fact that integrals of finite element 
functions must be evaluated by some numerical quadrature. In light of this, it should be noted that the 
finite dimensional space Vt on the evolving domain will not in general consist of piecewise polynomial 
functions. Furthermore, simplexes on the evolving domain will not in general be affine transformations of 
the reference simplex (see Figure [T] for an example). If formulation ( |3.18| ) is used and domain evolution 



is not spatially linear, the influence of numerical quadrature on the accuracy of the scheme should be 
considered. We leave this extension for future studies. 

7. Numerical experiments 

We now provide numerical evidence to back-up the estimate of Theorem[5]T] We use as a test problem, 
the Schnakenberg kinetics, although any other reaction kinetics that fulfils our assumptions could have 
been used. For the implementation we make use of the toolbox ALBERTA |23 | . The graphics were 
generated with PARAVIEW (10. 



7.1. Remark (Existence of solut ions to Problem [LI] with spatially linear isotropic evolution). In 

1301 , we showed that Problem [lT] equipped with the Schnakenberg reaction kinetics and Qt ^ C^i^t) is 
well posed under any bounded spatially linear isotropic evolution of the domain. If we assume this result 
holds on polygonal domains, we hav e sufficient regularity on the continuous problem to apply Theorem 



5.1 



and thus conclude scheme ( 3.16 ) with finite elements converges with optimal order. 

7.2. Definition (Experimental order of convergence). We denote the L^o (O, T; L2 {ft)^) error in the 
numerical scheme on a series of uniform refinements of a triangulation {c^} ._q ^ by {e^}i=o,...,Ar. 
The experimental order of convergence (EOC) is defined to be the numerical measure of rate of conver- 
gence of the scheme as /i^ ^ 0, where hn denotes the maximum mesh-size of ^ and is given by 

ln(ei+i/ei) 



(7.1) EOQ(e,,,+i,/i,,,+i) 



13. Numerical verification of th e a priori convergence rate. We examine the EOC of scheme ( |3.16| ) 



to approximate the solution to ( |l.l| with P^, and P^ basis functions and with uniform timestep r ~ /i^, 
T ^ It' and r ^ respectively (since the scheme is first order in time). 

We consider two different forms of domain evolution that include both domain growth and contraction 
to illustrate the versatility of the proposed finite element scheme. 

• Spatially linear periodic evolution: 

(7.2) A(0=C| 1 

• Spatially nonlinear periodic evolution: 

(7.3) {Atm^=^^ (l + ^sinf^V, 




T J 

for z = 1, . . . , (i. 

In both cases we take a time interval of [0, 1], the initial domain as the unit square and the parameter 
n = 1. Thus in both cases the domain grows to a square of length 2 at t = 0.5 before contracting back 
to the initial domain at end time. We take the diffusion coefficients D = (0.01, 1)^ and the parameter 
7 = 1. Problem [1 . 1 1 equipped with nonlinear reaction kinetics does not admit any closed form solutions. 
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In order to provide numerical verification of the convergence rate, we insert a source term such that the 
exact solution is. 



(7.4) 



U2 (C, t) = - sin(7rt) exp(-10 



Tables [T] and [2] show the EOCs for the two benchmark examples. For the first benchmark example 
where domain growth is linear with respect to space, we assemble the linear systems on the evolving 
domain corresponding to scheme ( |6.11 ). For the second benchmark example, as domain evolution is 
nonlinear, we assemble the linear systems on the reference frame corresponding to scheme (6.2). In both 
cases we observe that the error converges at the expected rate of ^ + 1 where £ is the degree of the finite 



element basis functions, providing numerical evidence for the a priori estimate of Theorem 5.1 





h 


2-5/2 


2-3 


2-7/2 


2-4 




e 


.119869 


.079248 


.040833 


.020521 




EOC 




1.1940 


1.9132 


1.9853 




e 


.035835 


.012380 


.004405 


.001567 




EOC 




3.0667 


2.9816 


2.9822 




h 


2-3/2 


2-2 


2-5/2 




p3 


e 


.153871 


.042748 


.010899 


.002728 




EOC 




3.6956 


3.9433 


3.9966 



Table 1 . Error in the L^o (O, T; L2 (^t)^) norm and EOCs for a benchmark problem 
with spatially linear domain evolution ( |7.2|). T he linear systems are assembled on the 
evolving frame corresponding to systeni^6JT]). 





h 


2-5/2 


2-3 


2-7/2 


2-4 




e 


.043839 


.022429 


.011133 


.005482 




EOC 




1.9337 


2.0210 


2.0441 


p2 


e 


.015906 


.005659 


.002005 


.000709 




EOC 




2.9822 


2.9941 


3.0010 




h 


2-3/2 


2-2 


2-5/2 


2-3 


p3 


e 


.022086 


.005514 


.001375 


.000344 




EOC 




4.0038 


4.0068 


4.0007 



Table 2. Error in the Lqo ^0, T; L2 norm and EOCs for a benchmark prob- 

lem with nonlinear domain evolution ( |7.3| ). The linear systems are assembled on the 
reference frame corresponding to system (|6.2|). 
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Reaction kinetic parameters Growth and numerical parameters 

Di D2 a b n T r DOFs 

.01 1.0 0.1 0.1 0.9 4 2000 10"^ 8321 
Table 3. Parameter values for numerical experiments with the Schnakenberg kinetics. 



7.4. Spatial pattern formation on evolving domains. We now present numerical results illustrating the 
influence of domain evolution on pattern formation by RDSs. We first present results for the Schnaken- 



berg kinetics with domain growth functions ( |7.2[ ) and ( |7.3| ) with identical initial conditions and identical 
numerical and reaction kinetic parameter values as given in Table [3] 

We take the unit square as the initial domain and the evolution of the boundary curve of the domain 
is identical in both cases, with the domain growing from a square of length 1 to a square of length 5 at 
t = 1000 before contracting to a square of length 1 at final time. In Figures [2] and [3] we show snapshots 
of the discrete activator (I^i) profiles obtained using the different schemes. The substrate profiles (W2) 
have been omitted as they are 180° out of phase with those of the activator. The snapshots in Figure 



[2] were computed using the evolving domain formulation ( 6.11| ), whilst the snapshots in Figure [3] were 



computed using the reference domain formulation ( |6.2[ ). Considering the Figures |2] and |3j we observe 
that no patterns are expressed when the domain is at initial or final times. This is as expected since linear 
stability analysis on a fixed square with the same length as the domain at t = and t = T indicates 
that, for the set of parameter values selected, the domain is below the critical domain size for which 
there exist admissible patterns. In the spatially linear case (Figure [2]), an initial half spot pattern forms 
which continuously transitions as the domain grows into a single spot positioned in the centre of the 
domain. As the domain contracts this single spot disappears (via spot annihilation) with the final domain 
exhibiting no spatial patterning. In the spatially nonlinear growth case (Figure [3]) the pattern transition 
is completely different with a half spot forming which splits as the domain grows to form two half spots 
which move around the domain before being annihilated as the domain contracts. The difference in 
patterning observed appears to be due to the differences in the growth function, with the results clearly 
illustrating the robustness of the numerical method in dealing with the complex non-uniform forms of 
domain evolution that are likely to be encountered in the biological problems we have in mind. 
We next consider the same Schnakenberg system on a domain with evolution of the form 

(7.5) {Atmi = 

with final time T = 1000. We select the same parameter values as Table |3] apart from the parameter 
7 which we set equal to 1. The initial domain (and reference domain) is taken as the square [—1,1]^. 
Figure]?] shows the activator concentrations on a domain with evolution of the form ( |7.5| ). We observe the 
now familiar spot- splitting behaviour as the domain grows and spot- annihilation as the domain contracts. 
Interestingly, we observe that the spots appear to orient themselves to maintain a relatively uniform level 
of separation (which is a characteristic of Turing patterns due to their intrinsic wavelength) even with this 
highly nonlinear form of domain evolution. 

Along with 2 component RDSs, 3 component RDSs where a second inhibitor quenches established 
maximums have been widely studied ll20l . Many interesting phenomena, which do not occur in 2 compo- 
nent systems, such as out of phase oscillations and space-time patterning which does not reach a steady 
state (even on fixed domains) are observed. The applications of such 3 component systems are of much 
importance, for example in the modelling of cell polarisation during chemotaxis |19|. To illustrate the 
versatility of our method in dealing with multiple component systems, we now present results for a three 
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Figure 2. Snapshots of the 
discrete activator (ui) profile 
for the Schnakenberg reaction 
kinetics on a domain with spa- 
tially linear evolution at times 
0, 580, 1000, 1500 and 2000. 
The substrate (^^2) profile is 
omitted as they are out of phase 
with those of the activator. For 
parameter values see Table |3] 
The results are computed us- 
ing the evolving domain for- 
mulation ( |6.11| ) with a uniform 
timestep and initial mesh. We 
observe the formation of a half 
spot which reorients to a sin- 
gle spot positioned in the cen- 
tre of the domain. As the do- 
main contracts this single spot 
is annihilated with the domain 
at end time exhibiting no pat- 
terns. 



Figure 3. Snapshots of the 
discrete activator (ui) profile 
for the Schnakenberg reaction 
kinetics on a domain with spa- 
tially nonlinear evolution at 
times 0, 590, 1000, 1750 and 
2000. For parameter values 
see Table |3l The results are 
computed using the reference 
domain formulation ( [6.2| ) with 
a uniform timestep and mesh. 
Although the evolution of the 
domain boundary is identical 
to the spatially linear case, the 
pattern transitions observed are 
markedly different with the for- 
mation of a half spot that splits 
into two half spots which are 
annihilated as the domain con- 
tracts to leave a final domain 
with no patterns. 
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(a) t = (b) t = 10 



(c) t = 30 



(d) t = 100 



(e) t = 190 





■M 



(h) t = 850 



(i) t = 950 



G) t = 1000 



Figure 4. Snapshots of the discrete activator (ui) profile for the Schnakenberg reac- 
tion kinetics on the reference domain and mapped to the evolving domain correspond- 



ing to system ( |6.2| ) and nonlinear domain evolution of the form ( |7.5| ). For parameter 
values see text. We see that spots form and increase in number (via splitting) as the 
domain grows and decrease in number (via annihilation) as the domain contracts. The 
spots appear to posses a common amplitude and maintain a relatively uniform separa- 
tion distance on the evolving domain. 



component model. We first consider a fixed domain, to illustrate both the qualitatively different solution 
behaviour to the 2 component case and that our method is readily applicable to multicomponent RDSs 
posed on fixed domains. 



7.5. Definition (Global and local inhibition kinetics [Wl). The following model is comprised of a sin- 
gle self-enhancing activator ui antagonised by global (fast diffusing) and local (slow diffusing) inhibitors; 
U2 and Us respectively: 



(7.6) 



s(n?+6i) 

7 {s' 



fs {ui,U2,us) = 7 {rsui - rsus) , 



r2U2) 
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Di D2 Ds J s ss ri r2 rs hi 
1.0 100 0.01 10"^ 0.005 0.8 0.005 0.008 0.001 0.01 
Table 4. Parameter values for numerical experiments with the 3 component kinetics. 



T T DOFs 
5 10-^ 8321 

Table 5. Numerical parameter values for experiments with the 3 component kinetics 
on a fixed domain. 



where < s, 61 , ri , S3 , r2 , ra , 7 < 00. We select the set of parameter values given in Tables |4]and[5] We 
take the domain to be the square = [—1,1]^. Snapshots of the discrete solutions are reported in Figure 
|5] An initial 4 spot pattern forms with a rapid transition to a 2 spot solution. The 2 spots then move around 
the domain boundary in a travelling wave like fashion. As expected from intuition, the activator ui and 
local inhibitor us concentrations are significantly more localised than the global inhibitor concentration 
Us. For details on the behaviour of 3 component systems (on fixed domains) and heuristic explanations 
of the observed behaviour we refer to |20|. The behaviour is markedly different to the 2-species case. No 
spatially inhomogeneous steady state is observed. None of the species are in or out of phase with spatial 
peaks of the global inhibitor travelling slightly ahead of spatial peaks of the activator which travels ahead 
of the spatial peaks of the local inhibitor. 

We now consider the same system posed on an evolving domain. Figure [6] shows snapshots of the 
evolution on a domain with growth of the form ( |7.2| ) and parameter values as given in Tables |4] and 
|6] Broadly speaking similar behaviour to the 2 component case is observed with the patterning mode 



T T DOFs 
4 100 10-^ 8321 
Table 6. Parameter values for numerical experiments with the 3 component kinetics 
on an evolving domain. 



(number of spots) generally increasing as the domain grows and decreasing as the domain contracts. The 
initial pattern is a 2 peak pattern that exhibits travelling wave behaviour similar to the fixed case. After 
the domain has grown sufficiently large new peaks appear either via splitting of existing peaks or peak 
insertion. When multiple peak patterns are present the behaviour is much more complicated than the 
two component case, the peaks travel around the domain and can collide leading to peak-merging or 
when peaks are in close proximity peak-annihilation, as reported in Figure |7] This behaviour of peak- 
merging and annihilation was observed previously by Venkataraman et al. 1 30 1 for 2 component systems, 
the novelty in the 3 component case is that this behaviour is observed even as the domain grows. This 
is of interest as with this specific 3 component system domain growth (contraction) does not lead to a 
monotonic increase (decrease) in the number of peaks. Our observations of these novel behaviours of 
3 component systems on evolving domains is to the best of our knowledge new and definitely warrants 
further investigation. 



Implicit-explicit timestepping with finite element approximation of reaction-diffusion systems on evolving domains 
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Figure 5. Snapshots of the discrete activator ui (left), discrete global inhibitor U2 
(middle) and discrete local inhibitor us (right) for the 3 component system ( |7.6| ) on 
a fixed domain at times 0.65, 1, 1.5, 2, 3.5 and 5. The colour (grey-scale) legend 
to the right of each column indicates the numerical range of each component during 
the experiment. The global inhibitor concentration is less spatially localised than the 
activator and local inhibitor concentration. We observe travelling wave like solutions. 
An initial 4 peak solution forms, subsequently 2 peaks are annihilated leaving a 2 peak 
solution which travels around the boundary of the domain for the remainder of the 
evolution. 
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Figure 6. Snapshots of the discrete activator ui (top), discrete global inhibitor U2 



(middle) and discrete local inhibitor us (bottom) for the 3 component system (7.6 ) on a 
domain with spatially linear evolution at times 0.5, 1, 2, 3, 4, 5, 6, 7, 8, 9 and 10 reading 
clockwise. The colour (grey-scale) legend above each row indicates the numerical 
range of each component during the experiment. Generally more peaks are observed 
on larger domains, however the behaviour is much richer than the 2 component case 
with the travelling wave like nature of the solutions leading to peak collisions which 
result in spot merging or peak annihilation when spots move into close proximity. 
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Figure 7 . Snapshots of the activator profiles of the 3 component system on the evolv- 
ing domain between t = 4.75 and t = 4.8 at intervals of 0.01. We observe peak 
merging, with the two peaks in the bottom left hand corner (within the dashed oval) 
merging and peak annihilation with the right most of the two peaks contained in the 
dashed rectangle being annihilated. Note the mode number of the pattern (number of 
spots) ait = 4.8 is less than ait = 4.75. 
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